%varying lamdac
clear all
clc
muc=0.033;
C=40;
Lambdac= 2:0.1:6 ;
%Lambdac=3;

Lambdav=8;
n=ceil(((Lambdav/muc) -1)/C);
% Pv=Gaussiandist(n);
Pv=decreasingdist(n);
Pc=Gaussiandist(n);
var_num = ((n+1)*(n+2))/2; 

lb=zeros(var_num,1);
ub=ones(var_num,1);
ub(var_num)=Inf;
A1=zeros(1,var_num);
for i=1:n
A1(1,i)=-Lambdav*Pv(1,i);
end
b1=ones(1,1)*muc*n*C - Lambdav;
A2=zeros(1,var_num);
A2(1,1)=Lambdav*Pv(1,1);
b2=ones(1,1)*muc;
Alin=[A1;A2];
blin=[b1;b2];

Aeq=zeros(n,var_num);
beq=ones(n,1);
begins=n;
ends=n;
for i=1:n
    begins=ends+1;
   ends=ends+i;
    for j=begins:ends
        Aeq(i,j)=1;
    end
   
end

b=zeros(n,1);
Waitings=zeros(length(Lambdac),1);
R=zeros(length(Lambdac),1);
exitflags=7*ones(length(Lambdac),1);
options = optimset('MaxFunEvals',40000);
for j=1:length(Lambdac)
for i=1:n
    b(i)= Lambdac(j)*Pc(i) ;
end
ub(var_num)=(Lambdav-Lambdac(j))/n;
cond = @(x)parameterfun(x,n,var_num,Pv,Lambdav,b);
obj=@(x)x(var_num)*(-1);
x0=zeros(var_num,1);
%x0 = ones(var_num-1,1);
%x0(var_num)= (Lambdav-Lambdac(j))/n;
x_best=zeros(var_num,1);
    for q=1:20
        [x,fval,exitflag,output,lambda] = fmincon(obj,x0,Alin,blin,Aeq,beq,lb,ub,cond,options);
        %x(var_num)
        if x(var_num)> x_best(var_num)
            x_best=x;
            exitflag_best=exitflag;
            lambda_best=lambda;
            %best_result=[x,fval,exitflag,output,lambda];
            x0=x;
        end
    end
exitflags(j)=exitflag_best;
R(j)= (1/x_best(var_num));
%[cout,ceqout] = parameterfun(x,n,Pv,Lambdav,b);
%Waitings(j)=max(1/(-cout(1:n)+0.1));
end


plot(Lambdac,R)
hold on
%smoothening
% for j=1:(length(Lambdac)-1)
%     while R(j)>R(j+1) 
%     R(j)=(R(j-1)+R(j+1))/2;
%     end
% end
% for j=25:29
%     R(j)=R(j-1)+0.2;
% end
% for j=20:26
%     R(j)=R(j-1)+0.21;
% end